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ABSTRACT 

Intrinsic differences in the processes involved in the electromagnetic and hadronic 
shower development in the atmosphere have been evidenced by means of a 
careful analysis of the event image given by the ARGO-YBJ detector. The 
images have been analyzed at different length scales and their multifractal na- 
ture has been studied. The use of the multiscale approach together with a 
properly designed and trained Artificial Neural Network, allowed us to obtain 
a good discrimination power. If confirmed by further studies on different event 
cathegories, this result would allow to nearly double the detector sensitivity to 
gamma ray sources. 



1 Introduction 



Gamma ray astronomy at energies around 100 GeT^-10 TeV is the main sci- 
entific goal of the ARGO-YBJ experiment ^1 121) , The detector, which is now 



being assembled in Tibet (China) at 4300 m a.s.l., is a full coverage Extensive 
Air Shower array consisting of a Resistive Plate Chamber (RPC) carpet of more 
than 6000 m 2 . It is logically divided into 154 units called clusters (7.64 x 5.72 
m 2 ), made by 12 RPCs (see Fig(TJ). Each RPC (1.26 x 2.85 m 2 ) is read out by 
10 pads (62 x 56 cm 2 ), which are further divided into 8 different strips (62 x 7 
cm 2 ), which provide the highest available space resolution. The signals coming 
from all the strips of a given pad are sent to the same channel of a multihit TDC. 
The whole system is designed in order to provide a single hit time resolution 
at the level of 1 ns, thus allowing a complete and detailed three-dimensional 
reconstruction of the shower front. The high altitude ( ~ 606 g/cm 2 ) and 
the full coverage ensure a very low primary energy threshold (E 7 wlOO GeV), 
while the detector time resolution gives a good pointing accuracy, thus allowing 
a high sensitivity to 7-ray sources. 

Gamma/hadron discrimination is a key issue in Very High Energy (VHE) 
gamma ray astronomy since it allows, together with a good angular resolution, 
the rejection of the huge background due to charged primary hadrons. The 
use of a full coverage detector with a high space granularity - like ARGO- 
YBJ - can give detailed images of the shower front. Intrinsic differences in the 
processes involved in the electromagnetic and hadronic shower development in 
the atmosphere can then be evidenced by means of a careful analysis of the 
event ^ El . Recently the use of multiscale behavior of event images has been 
showed to give good results in experiments exploiting the Imaging Atmospheric 
Cherenkov Techniques (IACT), together with the use of the so-called Hillas 
parameters ^ ^ . In these experiments the image is integrated over the entire 
shower development, while in the case of ARGO-YBJ a section of the shower 
is provided at a given (slanted) depth only, thus giving a potential lack of 
information. However the typical disuniformities present in hadronic events 
might be better evidenced in the case of ARGO-YBJ , these being partially 
masked in the case of IACT detectors because of the smearing effect due to the 
integration of the information along the shower development. 

In this work event images have been analyzed at different length scales and 
their multifractal nature has been studied. In particular the Discrete Wavelet 
Transforms have been applied since they allowed a differential approach to 
multifractality, that gave a higher discrimination power. 



Figure 1: General layout of the ARGO-YBJ detector. Each indicated box 
represents a 7.64 x 5.72 m 2 cluster of 12 RPCs (see text). 

2 The multiscale approach 

2.1 MultiFractal Analysis (MFA) 

We have considered the shower image seen by ARGO-YBJ as a lunction <f>(x) 
defined on a two-dimensional space and corresponding to the amplitudes given 
by the measured strip multiplicity. As a first step we can calculate the multi- 
fractal moment Z q (£) of order q, defined (see d) at the length scale £ as: 

Z,(*) = $>(2,*) 9 (1) 

{S} 

where 

p(2,i) = J_ / <p(x')dx" (2) 

gives the probability for a hit to be in the box Bi(x) of dimension I centered 
at x, being N tot — J 4>(x') dx' the total content of the image. By dividing the 
image, at different steps, into non overlapping pixels of size £, we can evaluate 
the box-amplitudes p(x, £). As pointed out in El El 01 j a t each order q, the MF 
moment is expected to have a power law dependence on £ in the high resolution 
limit, namely: Z q (£) - £ T ^ for I 1. 

By fitting the behavior of Z q {£) on £ (for each value of q) the MF scaling 
exponent r(q) can be extracted. The dependence of r(q) on q gives the main 
information on the MF properties of the image. 



Figure 2: Examples of Haar wavelet functions w(x, 2£) for 1=1,2 



2.2 The Discrete Wavelet Transform Analysis (DWTA) 

As shown in 3 ; the MF approach might not sufficiently characterize the image. 
An approach to multifractality based on the discrete wavelet transformations 
(DWT) is more appropriate. DWT can be seen as an expansion of <j>(x) on a 
discrete set of basis functions that are generated by scaling a so-called mother 
wavelet. For sake of simplicity, let us consider a one-dimensional case. If the 
Haar mother wavelet w(x,£) is chosen (see Fig|5Jl, the differences between the 
box amplitudes of two adjacent cells can be given by the proper convolution: 

p(x, I) ~ p(x + £,£) = — !— / <f>(x') w{x - x, 21) dx (3) 
Ntot J 

In this simple case, the mother wavelet can be indicated as (+ -). A differential 
approach to multifractality is then given by defining the DWT moment as: 

W q (l) = Y,\p(x,e)-p(x + £,£){« (4) 

{x} 

Also in this case the scaling properties of the image can be evidenced at high 
resolutions: W q {£) - £^ when l^>\. 

In the case of a two-dimensional distribution (i.e. an image) we need 
three base wavelets that can be indicated as: 

(--) c=) (-;) m 

(1) (2) (3) 

This led to write three different DWT moments Wq % \i), i = 1,2, 3, defined as: 

W ^W= E \p(x,y^)+p(x + £,yJ)-p(x,y + £J)- P (x + £,y + £,£)]" 
{x,y} 



% (2) W = E \p(x,y,e)-p(x + t,y,i)+ P (x,y + e,e)- P (x + e,y + e,e)\« (6) 

{z,y} 

These functions are expected to scale like Wq ~ £^*^ q \ when £ — * 1, with dif- 
ferent exponents (3 l {q), whose dependence on q gives the maximum information 
on the image properties. 

3 The Monte Carlo Samples 

We have generated ~2.8T0 5 gamma-initiated showers and ^2.6T0 5 proton- 
initiated ones making use of the CORSIKA code . The events have been 
taken within the energy range 30 GeV-^100 TeV with azimuth between and 15 
degrees and core at the detector center. The primary energy spectrum has been 
generated according to the measured power laws with spectral index 7=2.5 for 
gammas and 7=2.7 for hadrons. 

The detector response has been fully simulated by using ARGOG, a tool 
(based on the GEANT3 package El) developed within the ARGO-YBJ col- 
laboration. Since it is the hit multiplicity which is actually measured in EAS 
experiments, it is correct and appropriate to classify events following this vari- 
able and not in terms of the primary energy. Therefore we divided the simulated 
data sample into five multiplicity windows as reported in Tab^ where the av- 
erage primary energies are also shown. The 7 and proton energies are obviously 
different since they depend on the processes involved in the shower develop- 
ment in the atmosphere. In particular proton-induced showers results in lower 
multiplicity events with respect to the case of photon-induced ones with the 
same primary energy. Moreover in order to avoid distortion effects due to the 
finite size of the energy window in which the samples have been simulated, 
we have considered in the analysis only events which produced more than 50 
and less than 6000 hits. This corresponds to average primary photon energies 
between E 1 ~ 500 GeV and E y ~ 10 TeV (see TabQJ. 

4 Event analysis 

The multiresolution quantities defined in Sec|21have been used to analyze each 
event image. As reported in Sec^the ARGO-YBJ detector is made by a central 
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Figure 3: Dependences of the scaling exponents on the moment q for a sample 
of gamma and proton initiated- showers with energy from 3 TeV to 10 TeV . 
The error bars refer to the r.m.s values of the variuos distributions. 



carpet and a guard ring (see Fig^). In order to preserve the same symmetry 
at different lenght scales £, we decided to neglect, in this first analysis, the 
information coming from the external ring. We also decided to mask the central 
carpet with a square grid. In particular, since it is made of (120x 130)pads, the 
first and the last row of pads were not considered, while four empty columns 
of pads were added (two on the left and two on the right), thus obtaining 
a (128x128) pad mask. In order to limit statistical fluctuations of the hit 
multiplicity in the smallest pixels, the minimum pixel size considered in the 
analysis was set at (2x2)pads - about 1.4 m 2 - which will then corresponds to 
the maximum resolution, i.e. £—1. The analysis of an event goes then throught 
different steps, each correponding to different lenght scales. At the n-th step 
(with n — 1,2,..., 6), the image is divided into 2™x2 n square pixels of size 
&n = f» j containing each 4£ 2 pads, and the total strip multiplicity is computed 
in each considered pixel. 

The values of Log(Z g (£)) and Log(W 9 (£)) have been calculated following 
Eq^Jand Eq[f)| and their dependences on Log(£) have been fitted with a first 
order polynomial in the region where the scaling is expected (i.e. £ — > 1 ), 
for different values of the moment q. The scaling exponents r(q) and /3(q), 
defined in Sec|2 have then been obtained for each event. Their dependences 
on the moment q are shown in Fig|3]for a sample of gamma and proton initiated 
events. As can be seen there is a separation of the average values of the scaling 
exponents between electromagnetic and hadronic showers, which is however 



partly masked by the large fluctuations. Therefore the separations of the MF 
and DWT parameters are not sufficient in order to give a good discrimination 
between e.m. and hadronic showers, unless an Artificial Neural Network (ANN) 
is used as in Wc then decided to use multifractal parameters as inputs 

to a properly designed and trained ANN (see SecEJ. 

In order to increase the 7/h separation, a study on the shape and the 
simmetry of the event image has also been made. In particular we studied the 
skewness of each event by means of the third moment of the distributions of 
the hit coordinates on the detector plane, namely x and y. In general, the 
skewness of a sequence of values S = x%, X2, x%, Xn is defined as: 

where x and s are the average and r.m.s. over S, respectively. In our case we 
studied the behaviour of the quantity £ = x cu be/y C ube, with: 

_ E» "is? _ E a ntVi rS N 

Ei n * E* n i 

where n$ is the strip multiplicity of the pad at the position (xi, j/j). The value 
of £ is reported in Fig^for two samples of e.m. and hadronic showers in two 
different multiplicity windows. As can be seen, its average is almost one as 
expected from both its definition and the detector geometry, while the r.m.s is 
always larger for proton-initiated showers, due to the large fluctuations present 
in hadronic events. This behaviour suggested the use of £ as one of the inputs 
to the ANN. 



5 The Artificial Neural Network 



In order to perform the 7/hadron discrimination using the parameters we have 
introduced in the previous sections, we decided to make use of an Artificial 
Neural Network. The Neural Network we have chosen is of the feed forward type 
and it is made of 3 perceptrons layers. The ANN input is an eight-dimensional 
vector whose elements are: the event total hit multiplicity Nhu, the value of f, 
the multifractal exponents r and (3 for q = 4, 6, 8. The output vector is defined 
in a one dimensional space: it is trained to be 1 for gamma-initiated events 
and for hadronic showers. 
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Figure 4: Distributions of £ for simulated events in a couple of multiplicity 
regions (see text). 



Networks were implemented and optimized by using the Stuttgart Neural 
Network Simulator (SNNS) tool ^ ' . In designing the ANN, its characteristics 
were been deeply studied in order to reach a good compromise between the 
increase of recognition capability and the processing time. The network training 
was separately performed in 5 multiplicity windows by using several thousands 
events (see TabQ]). The ANN were then tested by using an independent reduced 
sample of events and the 7 recognition efficiency e 7 , together with the proton 
contamination (1 — e p ), were measured. The quantities e 7 and e p reached their 
plateau value after few thousands of ANN training epochs. An example of 
ANN output for a couple of multiplicity windows is shown in Fig|SJ 

The detector sensitivity to 7-ray sources is defined as S — N^/y/Nh, 
where _/V 7 is the number of gamma-initiated events, while is the hadron 
contamination of the considered sample. The use of a 7/h discrimination tool 
like the one we are considering in this work makes the sensitivity S to be 
multiplied by the factor Q = e 1 j\J\ — Eh- 

In this work values of Q ~2 have been reached (see TabQJ, which are 
among the largest obtained with various techniques in the experiments working 
in the field K>. This result would allow to nearly double the ARGO-YBJ 
sensitivity to a given source obtained from the pointing accuracy alone ^3 or, 
equivalently, to reduce by a factor four the time needed to observe it above the 
hadron background, with a given statistical significance. 




Network Output : 800<nhit<1500 




Figure 5: Outputs of the neural network in two of the five considered multiplicity 
regions. 



nhits 


events (7) 


events (p) 


<E P > 


<E 7 > 


Q 


50-^100 


6657 


3862 


0.8 TcV 


0.5 TcV 


1.28±0.01 


100-^500 


11556 


6862 


1.8 TcV 


1.1 TeV 


1.42±0.02 


500^800 


2571 


1644 


4.9 TcV 


2.9 TeV 


2.01±0.10 


800^1500 


3087 


1963 


7.6 TcV 


4.6 TeV 


1.78±0.07 


1500-^6000 


4329 


3053 


18.4 TcV 


11.3 TcV 


1.78±0.06 



Table 1: Main characteristics of the simulated data sample (no. of ANN train- 
ing events, average primary energy, ....) together with the values of Q forj/h 
discrimination that resulted from this work. 



6 Concluding remarks 

We presented the first results of a multiscale image analysis performed on Monte 
Carlo events by taking into account all the processes of shower development 
in the atmosphere and a full simulation of the ARGO-YBJ detector response. 
The images have been analyzed at different length scales and their multifractal 
nature has been studied. A set of eight image parameters has been identified 
and used as the input for an Artificial Neural Network, which was then trained 
in order to discriminate gamma initiated from proton-initiated showers. 

Since this is the first attempt of this kind of analysis in a EAS detector 
like ARGO-YBJ , we decided to restrict this study to events with the core at 



the detector center and azimuth angles not larger than 15 degrees, while all 
the energies with the correct spectral dependencies have been simulated in a 
wide range. At this level we also neglected the contribution given by primary 
nuclei heavier than protons. This is a good first order approximation because of 
the proton-dominated cosmic ray composition in the considered energy region. 
Furthermore heavier-nuclei-induced showers would produce event patterns with 
characteristics even more different from gamma-initiated ones. 

If the results obtained in this first study will be confirmed by a further 
analysis on the whole event categories (now in progress) , the detector sensitivity 
to a given source would nearly double or, equivalently, the time needed to 
observe it above the hadron background, with a given statistical significance, 
would be reduced by a factor four. 

The best performances in 7/h discrimination have been obtained for pho- 
ton primary energies in the few TeV range, while at higher energy this analysis 
might be well complemented by measuring the muon content of the shower 
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